JLA-ONA-2023

knitr::opts_chunk$set(echo = TRUE)

rm(list = ls())

#install.packages("ona", repos = c("https://cran.qe-libs.org", "https://cran.rstudio.org")) # only run once
#install.packages("tma", repos = c("https://cran.qe-libs.org", "https://cran.rstudio.org")) # only run once

most recent data Ayano shared

data <- read.csv("~/Rprojects/ayano/PES_010v6(selected).csv")

units <- c("Speaker", "Role", "TimePhase")

codes_1 <- c("Functions", "Aesthetics", "User", "Vision", "Prototyping") # design actions codes
codes_2 <- c("CA", "ALoK", "CSU", "GCA", "Projective","Regulative","Relational") # shared epistemic agency codes

hoo_rules <- conversation_rules(
  (TimePhase %in% UNIT$TimePhase & Role %in% UNIT$Role)
)

notes: accum_1 and set_1 are for design action codes; accum_2 and set_2 are for shared epistemic agency codes

accum_1 <-
  contexts(data, 
           units_by = units, 
           hoo_rules = hoo_rules) %>%
  accumulate_contexts(codes = codes_1, 
                      decay.function = decay(simple_window, window_size = 4),
                      return.ena.set = FALSE, norm.by = NULL)

Note: Not doing means rotation so I can compare everything in the same space

set_1 <-
  # model(accum_1,
  #       rotate.using = "mean",
  #       rotation.params = 
  #         list(engineers=accum_1$meta.data$Role=="Engineer",
  #              servicedesigners=accum_1$meta.data$Role=="ServiceDesigner"))
  
  model(accum_1)

set_1 <-
  # model(accum_1,
  #       rotate.using = "mean",
  #       rotation.params = 
  #         list(engineers=accum_1$meta.data$Role=="Engineer",
  #              servicedesigners=accum_1$meta.data$Role=="ServiceDesigner"))
  
  model(accum_1)

set_2 <-
  # model(accum_2,
  #       rotate.using = "mean",
  #       rotation.params = 
  #         list(engineers=accum_2$meta.data$Role=="Engineer",
  #              servicedesigners=accum_2$meta.data$Role=="ServiceDesigner"))
  
  model(accum_2)

global visual parameter (let’s make sure all the plots use the same level of multiplier)

set_2 <-
  # model(accum_2,
  #       rotate.using = "mean",
  #       rotation.params = 
  #         list(engineers=accum_2$meta.data$Role=="Engineer",
  #              servicedesigners=accum_2$meta.data$Role=="ServiceDesigner"))
  
  model(accum_2)

Model 1: design action codes

Plot 1a: Overall Group Comparisons

traces = c(2:10)
plot(set_1, title = "Groups") |>
  units(
    points= set_1$points$Role$ProductDesigner, 
    point_position_multiplier = point_position_multiplier,
    points_color = c("red"),
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
   units(
    points= set_1$points$Role$Engineer, 
    point_position_multiplier = point_position_multiplier,
    points_color = c("blue"),
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
   units(
    points= set_1$points$Role$ServiceDesigner, 
    point_position_multiplier = point_position_multiplier,
    points_color = c("green"),
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  nodes(node_size_multiplier = 0.3,
            node_position_multiplier = node_position_multiplier,
            self_connection_color = c("black"))|>
      plotly::layout(showlegend = TRUE, legend = list(x = 100, y = 0.9)) |>
      plotly::style(name = "Product Designer", traces = traces[1]) |>
      plotly::style(name = "Engineer", traces = traces[2]) |>
      plotly::style(name = "Service Designer", traces = traces[3])

NA

Plot 1b: Group subtractions

prod_pts = set_1$points$Role$ProductDesigner 
eng_pts = set_1$points$Role$Engineer
serv_pts = set_1$points$Role$ServiceDesigner

prod_mean_net = set_1$line.weights %>% dplyr::filter(Role == "ProductDesigner") %>% colMeans()
eng_mean_net = set_1$line.weights %>% dplyr::filter(Role == "Engineer") %>% colMeans()
serv_mean_net = set_1$line.weights %>% dplyr::filter(Role == "ServiceDesigner") %>% colMeans()

plot(set_1, title = "Product Designer vs Engineer") |>
  units(
    points = prod_pts, 
    points_color = "red",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points = eng_pts, 
    points_color = "blue",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = prod_mean_net - eng_mean_net, # optional multiplier to adjust for readability
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("red","blue")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("red","blue"))

plot(set_1, title = "Product Designer vs ServiceDesigner") |>
  units(
    points = prod_pts, 
    points_color = "red",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points =serv_pts, 
    points_color = "green",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = prod_mean_net - serv_mean_net, # optional multiplier to adjust for readability
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("red","green")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("red","green"))

plot(set_1, title = "Service Designer vs Engineer") |>
  units(
    points = serv_pts, 
    points_color = "green",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points = eng_pts, 
    points_color = "blue",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = serv_mean_net - eng_mean_net, # optional multiplier to adjust for readability
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("green","blue")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("green","blue"))

Plot 2: All Means by phase

#get list of time phases and groups loop over or use map function

phases = unique(set_1$points$TimePhase)

groups = unique(set_1$points$Role)

col_list_2 = c("red",#phase 1
              "blue", #phase 2
              "green") #phase 3

traces = c(2:10)

x = plot(set_1, title = "Group Means by Phase")
count = 1
for (i in 1:length(phases)){
  for (j in 1:length(groups)){
    points = set_1$points %>% filter(TimePhase == phases[i], Role == groups[j])
    x = x |>
      units(
        points = points,
        point_position_multiplier = point_position_multiplier,
        points_color = col_list_2[i],
        show_mean = TRUE, show_points = F, with_ci = FALSE
        ) |>
      nodes(node_size_multiplier = 0.3,
            node_position_multiplier = node_position_multiplier,
            self_connection_color = c("black"))|>
      plotly::layout(showlegend = TRUE, legend = list(x = 100, y = 0.9)) |>
      plotly::style(name = paste0(phases[i]," - ",groups[j]), traces = traces[count])
    count = count + 1
  }
}

x

NA

Statistical Tests:

In these data, the observations are nested within participants, so we need to check if this nesting has a significant effect. We do this by constructing a confidence interval around the intraclass correlation coefficient (https://en.wikipedia.org/wiki/Intraclass_correlation). If the interval contains zero, nesting is not significant.

reg_dat = set_1$points %>% filter(ENA_DIRECTION == "response")

ICC::ICCest(Speaker,SVD1,reg_dat) #not significant

Warning: 'x' has been coerced to a factor

$ICC
[1] 0.05242903

$LowerCI
[1] -0.2098135

$UpperCI
[1] 0.5007296

$N
[1] 10

$k
[1] 3.272727

$varw
[1] 0.3002277

$vara
[1] 0.01661158

ICC::ICCest(Speaker,SVD2,reg_dat) #not significant

Warning: 'x' has been coerced to a factor

$ICC
[1] -0.1818112

$LowerCI
[1] -0.3333495

$UpperCI
[1] 0.1964787

$N
[1] 10

$k
[1] 3.272727

$varw
[1] 0.07137526

$vara
[1] -0.01098045

Nesting is not significant so we proceed with an OLS regression. Below we set up two models for each dimension and test whether including interaction effects in the analysis significantly improves the modes.

mod_x = lm(SVD1 ~ as.factor(TimePhase) + Role, data = reg_dat)
mod_y = lm(SVD2 ~ as.factor(TimePhase) + Role, data = reg_dat)

mod_x_2 = lm(SVD1 ~ as.factor(TimePhase)*Role, data = reg_dat)
anova(mod_x,mod_x_2)

Analysis of Variance Table

Model 1: SVD1 ~ as.factor(TimePhase) + Role
Model 2: SVD1 ~ as.factor(TimePhase) * Role
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1     28 2.6341                              
2     24 1.7559  4   0.87815 3.0006 0.03851 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

mod_y_2 = lm(SVD2 ~ as.factor(TimePhase)*Role, data = reg_dat)
anova(mod_y,mod_y_2)

Analysis of Variance Table

Model 1: SVD2 ~ as.factor(TimePhase) + Role
Model 2: SVD2 ~ as.factor(TimePhase) * Role
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)    
1     28 1.1117                                 
2     24 0.2156  4   0.89614 24.939 3.02e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Interactions are significant, so we should use the interaction models. Let’s view those models

summary(mod_x_2)

Call:
lm(formula = SVD1 ~ as.factor(TimePhase) * Role, data = reg_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.85356 -0.01457  0.00790  0.08624  0.32816 

Coefficients:
                                          Estimate Std. Error t value Pr(>|t|)   
(Intercept)                                -0.1652     0.1352  -1.221  0.23385   
as.factor(TimePhase)2                       0.7013     0.2066   3.395  0.00239 **
as.factor(TimePhase)3                       0.6249     0.2066   3.025  0.00585 **
RoleProductDesigner                        -0.5044     0.1913  -2.637  0.01443 * 
RoleServiceDesigner                        -0.6041     0.1913  -3.158  0.00425 **
as.factor(TimePhase)2:RoleProductDesigner   0.1937     0.2815   0.688  0.49803   
as.factor(TimePhase)3:RoleProductDesigner   0.6236     0.2922   2.134  0.04324 * 
as.factor(TimePhase)2:RoleServiceDesigner  -0.2347     0.2815  -0.834  0.41272   
as.factor(TimePhase)3:RoleServiceDesigner   0.6444     0.2815   2.289  0.03117 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2705 on 24 degrees of freedom
Multiple R-squared:  0.8261,    Adjusted R-squared:  0.7681 
F-statistic: 14.25 on 8 and 24 DF,  p-value: 2.062e-07

summary(mod_y_2)

Call:
lm(formula = SVD2 ~ as.factor(TimePhase) * Role, data = reg_dat)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.196826 -0.035322 -0.001122  0.020930  0.273378 

Coefficients:
                                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                               -0.17794    0.04739  -3.755 0.000977 ***
as.factor(TimePhase)2                     -0.08890    0.07239  -1.228 0.231326    
as.factor(TimePhase)3                      0.62627    0.07239   8.651 7.68e-09 ***
RoleProductDesigner                        0.24264    0.06702   3.620 0.001366 ** 
RoleServiceDesigner                        0.29823    0.06702   4.450 0.000168 ***
as.factor(TimePhase)2:RoleProductDesigner  0.06701    0.09865   0.679 0.503477    
as.factor(TimePhase)3:RoleProductDesigner -0.35748    0.10237  -3.492 0.001880 ** 
as.factor(TimePhase)2:RoleServiceDesigner -0.27608    0.09865  -2.799 0.009962 ** 
as.factor(TimePhase)3:RoleServiceDesigner -0.93799    0.09865  -9.508 1.30e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09478 on 24 degrees of freedom
Multiple R-squared:   0.89, Adjusted R-squared:  0.8534 
F-statistic: 24.28 on 8 and 24 DF,  p-value: 1.032e-09

Viewing interactions

emmip(mod_x_2, Role ~ TimePhase)

emmip(mod_y_2, Role ~ TimePhase)

When interactions are introduced, interpreting regression coefficients becomes more difficult, so we calculated the marginal means and various contrasts (i.e., comparisons between groups of interest) to address our research questions.

RQ1 (X dimension - Design Actions): Here we want to compare the “main effect” of Role on the x dimension—i.e., are there differences between the roles on the X dimension when averaging over TimePhases.

emm_x1 = emmeans(mod_x_2, specs = pairwise ~ Role, weights = "proportional")

NOTE: Results may be misleading due to involvement in interactions

print(emm_x1$contrasts)

 contrast                          estimate    SE df t.ratio p.value
 Engineer - ProductDesigner           0.251 0.119 24   2.116  0.1077
 Engineer - ServiceDesigner           0.487 0.116 24   4.192  0.0009
 ProductDesigner - ServiceDesigner    0.236 0.113 24   2.085  0.1142

Results are averaged over the levels of: TimePhase 
P value adjustment: tukey method for comparing a family of 3 estimates 

Only significant difference on the X is between Engineer and Service designer. Let’s calculate the effect size (difference in standard deviations) of that difference. We will use Cohen’s d.

cohensd = function(diff_,g1,g2){
  diff_/(sqrt((sd(g1)^2 + sd(g2)^2)/2))
}

contrasts = as_tibble(emm_x1$contrasts)
diff_ = contrasts$estimate[2]
g1 = reg_dat %>% filter(Role == "Engineer") %>% select(SVD1)
g2 = reg_dat %>% filter(Role == "ServiceDesigner") %>% select(SVD1)

cohensd(diff_ = diff_,g1 = g1$SVD1,g2 = g2$SVD1)

[1] 0.9779281

RQ1 (Y dimension - Design Actions): Here we want to compare the “main effect” of Role on the y dimension—i.e., are there differences between the roles on the X dimension when averaging over TimePhases.

emm_y1 = emmeans(mod_y_2, specs = pairwise ~ Role, weights = "proportional")

NOTE: Results may be misleading due to involvement in interactions

print(emm_y1$contrasts)

 contrast                          estimate     SE df t.ratio p.value
 Engineer - ProductDesigner          -0.157 0.0415 24  -3.771  0.0026
 Engineer - ServiceDesigner           0.078 0.0407 24   1.917  0.1556
 ProductDesigner - ServiceDesigner    0.235 0.0397 24   5.915  <.0001

Results are averaged over the levels of: TimePhase 
P value adjustment: tukey method for comparing a family of 3 estimates 

Significant differences between engineer/product designer and productdesigner/servicedesigner. Let’s get the effect sizes:

contrasts = as_tibble(emm_y1$contrasts)

diff_evp = contrasts$estimate[1]
g1 = reg_dat %>% filter(Role == "Engineer") %>% select(SVD1)
g2 = reg_dat %>% filter(Role == "ProductDesigner") %>% select(SVD1)

d_evp = cohensd(diff_ = diff_,g1 = g1$SVD1,g2 = g2$SVD1)

diff_pvs = contrasts$estimate[3]
g1 = reg_dat %>% filter(Role == "ProductDesigner") %>% select(SVD1)
g2 = reg_dat %>% filter(Role == "ServiceDesigner") %>% select(SVD1)

d_pvs = cohensd(diff_ = diff_,g1 = g1$SVD1,g2 = g2$SVD1)

d_evp

[1] 0.9137988

d_pvs

[1] 0.8072687

Marginal means - x RQ2

Here we use marginal means to test for differences between timephases within groups and groups within timephase and

library(emmeans)

emm_x1a = emmeans(mod_x_2, specs = pairwise ~ TimePhase|Role, weights = "proportional") #timephases w/n groups
emm_x1a = as_tibble(emm_x1a$contrasts)

emm_x2 = emmeans(mod_x_2, specs = pairwise ~ Role|TimePhase, weights = "proportional") #groups w/n timephases
emm_x2 = as_tibble(emm_x2$contrasts)

#see here: https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html

emm_x1a = emm_x1a %>% filter(p.value < 0.05)
emm_x2 = emm_x2 %>% filter(p.value < 0.05)

emm_x1a

emm_x2

NA

Several contrasts are significant. Let’s calculate the effect sizes

print(list(t1_evp,t1_evs,t2_evs,t2_pvs))

[[1]]
[1] 1.826995

[[2]]
[1] 2.233997

[[3]]
[1] 4.14715

[[4]]
[1] 1.169235

Some of these are too large to be trusted. This is happening because some individuals are highly similar in the data. We will have to acknowledge this as it is likely due to small sample size

Marginal means -Y RQ2: Let’s repeate what we did for X here

emm_y1 = emmeans(mod_y_2, specs = pairwise ~ TimePhase|Role, weights = "proportional")
emm_y1 = as_tibble(emm_y1$contrasts)

emm_y2 = emmeans(mod_y_2, specs = pairwise ~ Role|TimePhase, weights = "proportional")
emm_y2 = as_tibble(emm_y2$contrasts)

emm_y1 = emm_y1 %>% filter(p.value < 0.05)
emm_y2 = emm_y2 %>% filter(p.value < 0.05)

emm_y1

emm_y2

Several contrasts are significant. Let’s calculate the effect sizes

print(list(t1_ep_y,t1_es_y,t2_ep_y,t2_ps_y,t3_es_y,t3_ps_y))

[[1]]
[1] -2.082427

[[2]]
[1] -2.734181

[[3]]
[1] -13.53816

[[4]]
[1] 2.197361

[[5]]
[1] 8.972244

[[6]]
[1] 8.239644

Model diagnostics

library(performance)
library(see)

check_model(mod_x_2)

Model has interaction terms. VIFs might be inflated.
  You may check multicollinearity among predictors of a model without interaction terms.

check_model(mod_y_2)

Model has interaction terms. VIFs might be inflated.
  You may check multicollinearity among predictors of a model without interaction terms.

Things look ok for the most part

Plot network subtractions (replace TimePhase and Role for different combinations)

mean1 = set_1$points %>% dplyr::filter(TimePhase == 1, Role == "Engineer") #update here
mean2 = set_1$points %>% dplyr::filter(TimePhase == 1, Role == "ServiceDesigner") #update here

mean1_net = set_1$line.weights %>% dplyr::filter(TimePhase == 1, Role == "Engineer") %>% colMeans()
mean2_net = set_1$line.weights %>% dplyr::filter(TimePhase == 1, Role == "ServiceDesigner") %>% colMeans() #update here

plot(set_1, title = "Network subtraction") |>
  units(
    points = mean1, 
    points_color = "red",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points = mean2, 
    points_color = "blue",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = mean1_net - mean2_net, # optional multiplier to adjust for readability
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("red","blue")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("red","blue"))

Model 2: Epistemic Agency Codes

Plot 1: Mean Network

node_size_multiplier = 0.3 # scale up or down node sizes
node_position_multiplier = 1 # zoom in or out node positions
point_position_multiplier = 1 # zoom in or out the point positions
edge_arrow_saturation_multiplier = 1.5 # adjust the chevron color lighter or darker
edge_size_multiplier = 2 # scale up or down edge sizes

plot(set_2, title = "Overall Mean--Agency") |>
  units(
    points= set_2$points, 
    point_position_multiplier = point_position_multiplier,
    points_color = c("black"),
    show_mean = TRUE, show_points = TRUE, with_ci = TRUE) |>
  edges(
    weights =set_2$line.weights,
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("black")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    node_labels = TRUE, # change this to FALSE can remove node labels in case you want to add them back in a nicer font or size for your presentations or publications
    self_connection_color = c("black"))

Plot 1b: Overall Group Comparisons

plot(set_2, title = "Groups") |>
  units(
    points= set_2$points$Role$ProductDesigner, 
    point_position_multiplier = point_position_multiplier,
    points_color = c("red"),
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
   units(
    points= set_2$points$Role$Engineer, 
    point_position_multiplier = point_position_multiplier,
    points_color = c("blue"),
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
   units(
    points= set_2$points$Role$ServiceDesigner, 
    point_position_multiplier = point_position_multiplier,
    points_color = c("green"),
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  nodes(node_size_multiplier = 0.3,
            node_position_multiplier = node_position_multiplier,
            self_connection_color = c("black"))|>
      plotly::layout(showlegend = TRUE, legend = list(x = 100, y = 0.9)) |>
      plotly::style(name = "Product Designer", traces = traces[1]) |>
      plotly::style(name = "Engineer", traces = traces[2]) |>
      plotly::style(name = "Service Designer", traces = traces[3])

NA

Plot 1c: Group subtractions

prod_pts = set_2$points$Role$ProductDesigner 
eng_pts = set_2$points$Role$Engineer
serv_pts = set_2$points$Role$ServiceDesigner

prod_mean_net = set_2$line.weights %>% dplyr::filter(Role == "ProductDesigner") %>% colMeans()
eng_mean_net = set_2$line.weights %>% dplyr::filter(Role == "Engineer") %>% colMeans()
serv_mean_net = set_2$line.weights %>% dplyr::filter(Role == "ServiceDesigner") %>% colMeans()

plot(set_2, title = "Product Designer vs Engineer") |>
  units(
    points = prod_pts, 
    points_color = "red",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points = eng_pts, 
    points_color = "blue",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = prod_mean_net - eng_mean_net, # optional multiplier to adjust for readability
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("red","blue")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("red","blue"))

plot(set_2, title = "Product Designer vs ServiceDesigner") |>
  units(
    points = prod_pts, 
    points_color = "red",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points =serv_pts, 
    points_color = "green",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = prod_mean_net - serv_mean_net, # optional multiplier to adjust for readability
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("red","green")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("red","green"))

plot(set_2, title = "Service Designer vs Engineer") |>
  units(
    points = serv_pts, 
    points_color = "green",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points = eng_pts, 
    points_color = "blue",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = serv_mean_net - eng_mean_net, # optional multiplier to adjust for readability
    edge_size_multiplier = 1,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("green","blue")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("green","blue"))

Plot 2: All Means

#get list of time phases and groups loop over or use map function

phases = unique(set_2$points$TimePhase)

groups = unique(set_2$points$Role)


col_list_2 = c("red",
              "blue",
              "green")

traces = c(2:10)

x_agency = plot(set_2, title = "Group Means by Phase -- Agency")
count = 1

for (i in 1:length(phases)){
  for (j in 1:length(groups)){
    points = set_2$points %>% filter(TimePhase == phases[i], Role == groups[j])
    #print(head(points))
    x_agency = x_agency |>
      units(
        points = points,
        point_position_multiplier = point_position_multiplier,
        points_color = col_list_2[i],
        show_mean = TRUE, show_points = F, with_ci = FALSE
        ) |>
      nodes(node_size_multiplier = 0.3,
            node_position_multiplier = node_position_multiplier,
            self_connection_color = c("black"))|>
      plotly::layout(showlegend = TRUE, legend = list(x = 100, y = 0.9)) |>
      plotly::style(name = paste0(phases[i]," - ",groups[j]), traces = traces[count])
    count = count + 1
    
  }
}

x_agency

Check nesting

reg_dat_agency = set_2$points %>% filter(ENA_DIRECTION == "response")

ICC::ICCest(Speaker,SVD1,reg_dat_agency) #not significant

Warning: 'x' has been coerced to a factor

$ICC
[1] 0.02333382

$LowerCI
[1] -0.2268987

$UpperCI
[1] 0.4706009

$N
[1] 10

$k
[1] 3.272727

$varw
[1] 0.2901337

$vara
[1] 0.00693167

ICC::ICCest(Speaker,SVD2,reg_dat_agency) #significant we will need to include a variable for individual in the model

Warning: 'x' has been coerced to a factor

$ICC
[1] 0.5156316

$LowerCI
[1] 0.1639302

$UpperCI
[1] 0.8233863

$N
[1] 10

$k
[1] 3.272727

$varw
[1] 0.08844428

$vara
[1] 0.09415284

mod_x_agency = lm(SVD1 ~ as.factor(TimePhase) + Role, data = reg_dat_agency)
mod_y_agency = lm(SVD2 ~ as.factor(TimePhase) + Role, data = reg_dat_agency)


summary(mod_x_agency)

Call:
lm(formula = SVD1 ~ as.factor(TimePhase) + Role, data = reg_dat_agency)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.72619 -0.10243  0.00017  0.18426  0.34417 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            -0.3092     0.1103  -2.804  0.00906 ** 
as.factor(TimePhase)2   0.6862     0.1201   5.713 3.97e-06 ***
as.factor(TimePhase)3   1.0597     0.1232   8.601 2.40e-09 ***
RoleProductDesigner    -0.2519     0.1257  -2.004  0.05486 .  
RoleServiceDesigner    -0.4308     0.1232  -3.496  0.00159 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2873 on 28 degrees of freedom
Multiple R-squared:  0.7564,    Adjusted R-squared:  0.7216 
F-statistic: 21.74 on 4 and 28 DF,  p-value: 2.999e-08

summary(mod_y_agency)

Call:
lm(formula = SVD2 ~ as.factor(TimePhase) + Role, data = reg_dat_agency)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.4179 -0.2658 -0.1107  0.2753  0.5530 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)   
(Intercept)            0.12309    0.12362   0.996  0.32793   
as.factor(TimePhase)2  0.01307    0.13464   0.097  0.92338   
as.factor(TimePhase)3 -0.05497    0.13812  -0.398  0.69364   
RoleProductDesigner   -0.50503    0.14093  -3.584  0.00127 **
RoleServiceDesigner    0.15828    0.13812   1.146  0.26151   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3221 on 28 degrees of freedom
Multiple R-squared:  0.4817,    Adjusted R-squared:  0.4076 
F-statistic: 6.505 on 4 and 28 DF,  p-value: 0.0007821

check for interactions: (interactions are significant so we should use those models)

mod_x_agency_2 = lm(SVD1 ~ as.factor(TimePhase) + Role + as.factor(TimePhase)*Role, data = reg_dat_agency)
mod_y_agency_2 = lm(SVD2 ~ as.factor(TimePhase) + Role + as.factor(TimePhase)*Role, data = reg_dat_agency)

summary(mod_x_agency_2)

Call:
lm(formula = SVD1 ~ as.factor(TimePhase) + Role + as.factor(TimePhase) * 
    Role, data = reg_dat_agency)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.72200 -0.04168  0.00319  0.12198  0.33028 

Coefficients:
                                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                               -0.46990    0.12236  -3.840 0.000788 ***
as.factor(TimePhase)2                      1.18141    0.18691   6.321 1.56e-06 ***
as.factor(TimePhase)3                      1.10001    0.18691   5.885 4.52e-06 ***
RoleProductDesigner                       -0.03821    0.17304  -0.221 0.827107    
RoleServiceDesigner                       -0.16249    0.17304  -0.939 0.357071    
as.factor(TimePhase)2:RoleProductDesigner -0.55246    0.25471  -2.169 0.040216 *  
as.factor(TimePhase)3:RoleProductDesigner -0.15842    0.26433  -0.599 0.554563    
as.factor(TimePhase)2:RoleServiceDesigner -0.84961    0.25471  -3.336 0.002761 ** 
as.factor(TimePhase)3:RoleServiceDesigner -0.00878    0.25471  -0.034 0.972787    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2447 on 24 degrees of freedom
Multiple R-squared:  0.8485,    Adjusted R-squared:  0.798 
F-statistic:  16.8 on 8 and 24 DF,  p-value: 4.228e-08

summary(mod_y_agency_2)

Call:
lm(formula = SVD2 ~ as.factor(TimePhase) + Role + as.factor(TimePhase) * 
    Role, data = reg_dat_agency)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.39436 -0.08520  0.00161  0.07954  0.63608 

Coefficients:
                                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                               -0.09125    0.11044  -0.826 0.416831    
as.factor(TimePhase)2                      0.58628    0.16871   3.475 0.001959 ** 
as.factor(TimePhase)3                      0.08627    0.16871   0.511 0.613779    
RoleProductDesigner                        0.03593    0.15619   0.230 0.820034    
RoleServiceDesigner                        0.26034    0.15619   1.667 0.108555    
as.factor(TimePhase)2:RoleProductDesigner -1.01499    0.22991  -4.415 0.000184 ***
as.factor(TimePhase)3:RoleProductDesigner -0.74981    0.23859  -3.143 0.004411 ** 
as.factor(TimePhase)2:RoleServiceDesigner -0.61494    0.22991  -2.675 0.013251 *  
as.factor(TimePhase)3:RoleServiceDesigner  0.23732    0.22991   1.032 0.312246    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2209 on 24 degrees of freedom
Multiple R-squared:  0.791, Adjusted R-squared:  0.7214 
F-statistic: 11.36 on 8 and 24 DF,  p-value: 1.662e-06

anova(mod_x_agency,mod_x_agency_2)

Analysis of Variance Table

Model 1: SVD1 ~ as.factor(TimePhase) + Role
Model 2: SVD1 ~ as.factor(TimePhase) + Role + as.factor(TimePhase) * Role
  Res.Df    RSS Df Sum of Sq     F  Pr(>F)  
1     28 2.3112                             
2     24 1.4373  4   0.87388 3.648 0.01853 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(mod_y_agency,mod_y_agency_2)

Analysis of Variance Table

Model 1: SVD2 ~ as.factor(TimePhase) + Role
Model 2: SVD2 ~ as.factor(TimePhase) + Role + as.factor(TimePhase) * Role
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1     28 2.9044                                  
2     24 1.1710  4    1.7333 8.8813 0.0001506 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Need to run mixed model for Y ##model is singular so stick with OLS

library(lmerTest)

#mod_y_agency_mixed = lmerTest::lmer(SVD2 ~ 1 + (1|Speaker) + as.factor(TimePhase)*Role ,data = reg_dat_agency)

#summary(mod_y_agency_mixed)

Marginal means x - RQ1

emm_ag1 = emmeans(mod_x_agency_2, specs = pairwise ~ Role, weights = "proportional")

NOTE: Results may be misleading due to involvement in interactions

print(emm_ag1$emmeans)

 Role             emmean     SE df lower.CL upper.CL
 Engineer         0.2572 0.0777 24    0.097   0.4175
 ProductDesigner -0.0131 0.0740 24   -0.166   0.1396
 ServiceDesigner -0.1911 0.0708 24   -0.337  -0.0449

Results are averaged over the levels of: TimePhase 
Confidence level used: 0.95 

print(emm_ag1$contrasts)

 contrast                          estimate    SE df t.ratio p.value
 Engineer - ProductDesigner           0.270 0.107 24   2.520  0.0476
 Engineer - ServiceDesigner           0.448 0.105 24   4.265  0.0008
 ProductDesigner - ServiceDesigner    0.178 0.102 24   1.737  0.2123

Results are averaged over the levels of: TimePhase 
P value adjustment: tukey method for comparing a family of 3 estimates 

Viewing interactions

emmip(mod_x_agency_2, Role ~ TimePhase)

Marginal means - x RQ2

emm_ag1 = emmeans(mod_x_agency_2, specs = pairwise ~ TimePhase|Role, weights = "proportional")
emm_ag1 = as_tibble(emm_ag1$contrasts)

emm_ag2 = emmeans(mod_x_agency_2, specs = pairwise ~ Role|TimePhase, weights = "proportional")
emm_ag2 = as_tibble(emm_ag2$contrasts)

emm_ag1 = emm_ag1 %>% filter(p.value < 0.05)
emm_ag2 = emm_ag2 %>% filter(p.value < 0.05)
                             
emm_ag1

emm_ag2

NA

Effect sizes

#todo

Marginal means y - RQ1

emm_ag_y1 = emmeans(mod_y_agency_2, specs = pairwise ~ Role, weights = "proportional")

NOTE: Results may be misleading due to involvement in interactions

print(emm_ag_y1$contrasts)

 contrast                          estimate     SE df t.ratio p.value
 Engineer - ProductDesigner           0.530 0.0968 24   5.470  <.0001
 Engineer - ServiceDesigner          -0.127 0.0949 24  -1.341  0.3868
 ProductDesigner - ServiceDesigner   -0.657 0.0925 24  -7.104  <.0001

Results are averaged over the levels of: TimePhase 
P value adjustment: tukey method for comparing a family of 3 estimates 

Viewing interactions

emmip(mod_y_agency_2, Role ~ TimePhase)

Marginal means - y RQ2

emm_ag_y1 = emmeans(mod_y_agency_2, specs = pairwise ~ TimePhase|Role, weights = "proportional")
emm_ag_y1 = as_tibble(emm_ag_y1$contrasts)

emm_ag_y2 = emmeans(mod_y_agency_2, specs = pairwise ~ Role|TimePhase, weights = "proportional")
emm_ag_y2 = as_tibble(emm_ag_y2$contrasts)

emm_ag_y1 = emm_ag_y1 %>% filter(p.value < 0.05)
emm_ag_y2 = emm_ag_y2 %>% filter(p.value < 0.05)

emm_ag_y1

emm_ag_y2

Calculate eff sizes

#to do

Model diagnostics

check_model(mod_x_agency_2)

Model has interaction terms. VIFs might be inflated.
  You may check multicollinearity among predictors of a model without interaction terms.

check_model(mod_y_agency_2)

Model has interaction terms. VIFs might be inflated.
  You may check multicollinearity among predictors of a model without interaction terms.

Looks ok

Plot network subtractions

mean1 = set_2$points %>% dplyr::filter(TimePhase == 1, Role == "Engineer")
mean2 = set_2$points %>% dplyr::filter(TimePhase == 1, Role == "ProductDesigner")

mean1_net = set_2$line.weights %>% dplyr::filter(TimePhase == 1, Role == "Engineer") %>% colMeans()
mean2_net = set_2$line.weights %>% dplyr::filter(TimePhase == 1, Role == "ProductDesigner") %>% colMeans()

plot(set_2, title = "Network subtraction") |>
  units(
    points = mean1, 
    points_color = "red",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points = mean2, 
    points_color = "blue",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = mean1_net - mean2_net, # optional multiplier to adjust for readability
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("red","blue")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("red","blue"))